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ABSTRACT 

The 21-cm hyperfine line of atomic hydrogen (Hi) is a promising probe of the cosmic 
dark ages. In past treatments of 21-cm radiation it was assumed the hyperfine level 
populations of H I could be characterized by a velocity-independent "spin temperature" 
T s determined by a competition between 21-cm radiative transitions, spin-changing 
collisions, and (at lower redshifts) Lya scattering. However we show here that, if the 
collisional time is comparable to the radiative time, the spin temperature will depend 
on atomic velocity, T s = T s (v), and one must replace the usual hyperfine level rate 
equations with a Boltzmann equation describing the spin and velocity dependence 
of the H I distribution function. We construct here the Boltzmann equation relevant 
to the cosmic dark ages and solve it using a basis-function method. Accounting for 
the actual spin-resolved atomic velocity distribution results in up to a ~2 per cent 
suppression of the 21-cm emissivity, and a redshift and angular-projection dependent 
suppression or enhancement of the linear power spectrum of 21-cm fluctuations of up 
to ~5 per cent. The effect on the 21-cm line profile is more dramatic — its full- width 
at half maximum (FWHM) can be enhanced by up to ~ 60 per cent relative to the 
velocity-independent calculation. We discuss the implications for 21-cm tomography 
of the dark ages. 

Key words: cosmology: theory - intergalactic medium - atomic processes - line: 
profiles. 



1 INTRODUCTION 

While both earlier and later epochs of the Universe have now been observed by astronomers, the cosmic dark ages (the time 
between the formation of the first atoms at z ~ 1000 and the subsequent formation of the first luminous objects) remain an 
observational enigma. The difficulty of probing this epoch stems from the relatively inert physical properties expected of the 
dark-age gas. It is diffuse, without the luminous galaxies or quasars that have proven so useful for lower-redshift astrophysics. 
It is mostly neutral, hence signatures such as recombination lines, free-free radiation, and Compton scattering of the cosmic 
microwave background (CMB) are essentially absent. Finally, its composition consists mainly of cold H and He atoms, which 
have very few lines that can be excited at the relevant temperatures. 

A notable exception to this last property is the hyperfine 21-cm line of neutral hydrogen (Hi). The Hi 21-cm line has 
attracted considerable attention as a possible probe of cosmic evolution during the dark ages, despite substantial observational 
challenges. In particular by using frequency as well as angular information, it should be possi ble to map out the 21-cm 
signal in three dimensions, thereby grea tly increasing the number of modes that can be observed jLoeb fc ZaldarriaealEool : 
Moral es et al j200StlBowman et al-bool) . Furthermore, the foregrounds to this signal are expected to have a distinctly different 
character in both angular and especially frequ ency space (they are expected to be relatively smooth in frequ ency space) which 
should allow for effective foreground removal ijZaldarriaga et al1l2004l : ISantos et alJl200nl : IWang et, alJl200Rl ) . 

The nature of the 21-cm signal depends on both the properties of the primordial gas and the mechanisms that excite 
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or de-excite the hydrogen atoms. The major mechanisms operating during the dark ages are the radiative transitions at 
A21 = 21.1cm and atom-atom collisions; the former tend to cause the hyperfine level populations to relax to the CMB 
temperature T 7 , whereas the latter tend to drive the level populations towards the gas kinetic temperature Tfc. The atomic 
gas is expected to be colder than the CMB during the dark ages since it cools adiabatically as T oc a" 2 after residual 
Compton heating becomes ineffe ctive. During this period w e thus expect to see the 21-cm signal in absorption when the 
collisional coupling is strong and lLoeb fc Zaldarriaeal i2004l) have shown that at high redshifts z ^ 50 the entire Universe 
should be collisionally coupled and an absorption signal should be present. At lower redshifts collisions become less efficient 
due to the low density and the level populations approach T 7 ; in this case the absorption signal fades because the Hi spins 
are nearly in thermal equilibrium with the CMB. Only the highest-density regions remai n collisionally coupled, and these 
may even appear in emission if shocks or adiabatic compression heat them to Tj, > T 7 Jlliev et al.ll2002l 120031) . At later 
times when sources of ultraviolet radiation have turned on, scattering of photons in the H I Lya resonance can also become 
significant, and can once again drive the level populations towards T/j even in mean-density and underdense regions. These 
early sources of rad iation may also have emitted X-rays, which could heat the intergalactic m edium (IGM) to well above the 
CMB temperature l|Madau et al.ll997tlchen fc Miralda-Escude1l2004l : lKuhlen fc Madaulboosh . and perhaps pro duce sufficient 
ioniz ation that atom-electron collisions are important in determining the spin temperature in some regions ((Kuhlen et, alJ 
l200rt . 

A detailed understanding of all these excitation and de-excitation mechanisms and their effect on the 21-cm emissivity 
and line profile will be necessary in order to interpret future 21-cm observations. In the 195 0s and 1960s, analysis of 21-cm 
radia tion from the in terstellar medium motivated pioneering s tudies of the Lya scattering JWouthuvsenlll952l : lFieldlll95ci 
1959) and collisional iDalgarnolll96ll: lAllison fc Dalgarnolll96il mechanisms. The results of these analyses have largely been 
carried over into the theory of 21-cm radiation from the high-redshift IGM and dark ages, with some improvements and 
adaptations to accommodate different physical situations. For example, several authors have computed the efficiency of the 
Lya scattering mechanism with detailed consideration of cosmological radiative transfer eff ects, interaction with Ly/3 and 
higher- order Hi Lyman lines, and the fine and hyperfine structure of the Lya reson ance llChen fc Miralda-Escudll2004l : 
IffiratalEoOfil: IPritchard fc Furlanettoll20o4 Ichuzhov fc Shapirolboosh . IZveelmanl I^OOSh has also re-evaluated the collisional 
hyperfine excitation rate using modern atomic physics techniques. 

In this paper we note that heretofore all studies of high-redshift 21-cm radiation have assumed that the hyperfine 
level population ratios of the hydrogen atoms in the primordial gas are independent of velocity. If the level populations 
are also isotropic, they dependent only on the total spin F and not on the magnetic quantum number Mp, so they can 
be characterized by a single spin temperature T s } This assumption implies that the ratio of hyperfine level populations is 
nm(F = l)/nui(F = 0) = 3e~ T *' /Ts independent of the velocity of the hydrogen atoms. 2 This should be the case if the level 
populations are determined solely by radiative transitions with the CMB in the 21-cm line or solely by atomic collisions as, in 
these cases, the hyperfine level populations will thermalize at a temperature T 7 or Tk respectively. However we demonstrate 
here that, because the level populations are determined by both radiative and collisional transitions and the collisional spin- 
change cross section for H-H collisions is velocity dependent while the radiative cross section is velocity independent, the 
hyperfine level-resolved distribution function of H I is not thermal. This is our main result and we show below that this has a 
significant (order unity) effect on the 21-cm line profile of collisionally coupled gas at high redshift. 

The paper is organized as follows: We outline the basic kinetic theory and derive the collisional piece of the linear 
Boltzmann equation in Section[5]and define the basis-function method we use to describe deviations from a thermal distribution 
function in Section [3] In Section [I] we evaluate the H-H collision term and in Section [K] we evaluate the He-H collision term 
to account for the effects of helium. We solve the Boltzmann equation in Section and present our results for the 21-cm line 
profile in Section^ A short discussion and summary of our main results follows in Section|H] We have included in AppendixlAl 
a summary of some useful properties of the Gauss- Hermite basis functions, and in Appendix [5] the details of the quantum 
mechanical calculation of the relevant atomic cross sections. When needed for our numerical results we have assumed a flat 
cosmology with matter density f2 m = 0.30, baryon density fib = 0.042, Hubble constant Ho = 72kms~ 1 Mpc -1 , primordial 
helium fraction Y Ke = 0.24, and CMB temperature T = 2.728 K. 



2 KINETIC THEORY 

In the standard calculation of the 21-cm signal one uses the populations nm(F — 0) and nm(F — 1) of the ground and 
hyperfine-excited levels as the basic quantities which describe the state of the hydrogen atoms. In this paper we aim to go 

1 Several g roups have discu ssed the possibility of Mp -dependent pop ulations due to magn etic fields JCoorav fc FurlanettclfcoOEft . CMB 
anisotropy JZvge lman 200^) , and a nisotropies in the Lya radiation iBabich fc Loebll2005l) . In all cases the effect was argued to be very 
small, and indeed Izygelmanl l2005l> showed that if any net spin polarization of the Hi atoms was initially present, it would be destroyed 
by radiative transitions much faster than the Hubble time. 

2 Here T* = hcj\i\kg = 68.2 mK is the energy of the 21-cm transition in temperature units. 
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beyond this calculation by solving for the joint distribution of spins and velocities. This will be represented by the distribution 
function, /f(i>), which describes the number of atoms in hyperfine level F per unit volume in position space, per unit volume 
in velocity space. For comparison, the standard calculation assumes that the distribution function can be separated into its 
kinetic degrees of freedom with temperature Tk and its spin degrees of freedom with temperature T 3 . In this case one can 
write 

/f(w)|t, = nmyF$T k (v), (1) 
where $T fc (u) represents the Maxwellian distribution, 

^W=(2^b-^, (2) 



<y = \/kBTk/rriH is the thermal velocity dispersion of the atoms, and v is the velocity of the atom in the reference frame 
of the baryonic fluid. Note that here and elsewhere in the paper the notation "|t s " is shorthand for "evaluated in the spin- 
temperature approximation." The factor uf in Eq. Q is the fraction of the atoms in the F hyperfine level, and an elementary 
calculation shows that this quantity is related to the spin temperature by 

- 1 1 3T» _ 3e- T * /T ° 3 _ 3?; 

V0 ~ 3e- r */ T =+l ~4 + 167; Vl ~ 3e- T */ T =+l ~ 4 " WF a ' ( ) 

where the approximation T* <C T s is valid in all cases of astrophysical interest. The standard analysis solves for T s using 
the continuity equation for F — and F = 1 Hi atoms. By determining all rates (collisional, radiative, and Lya) of 
production or destruction of Hi atoms in a given level one obtains the abundances n(F) and hence the spin temperature 
e ~ T */ T s — n(l)/[3n(0)]. In the standard case, the brightness temperature due to the 21-cm line is 



Tb\ Ts = 



H(z) dv\^ 

— lj. _| y — 

1 + z dr» 



3c 3 AT* ( T 7 ■ 

-F7T— m nHI ( 1 - W- ) » ( 4 ) 



32^,(1 + ^)2 ni V T, 



where H(z) is the Hubble rate, Wy P ° c ' is the radial component of the peculiar velocity, c is the speed of light, A is the 
21-cm Einstein coeffic ient, viq is the frequency of the 21-cm line, and T 7 is the CMB temperature feharadwai fc Alilbooi 
iBarkana fe Loebll2005h . One of the major goals of this paper is to update Eq. gj to include the full spin-velocity distribution. 



2.1 Simplifying approximations; distribution function 

The distribution function Jf(v) is a much more general description of the state of the gas than the three numbers (ni, Tk, T a } 
in the standard treatment. Nevertheless, the most general description possible would also account for dependence on Mf and 
quantum correlations among the various states of H I. In this section we describe the key simplifying approximations that 
allow us to fully characterize the Hi atoms by a distribution function. 

Since the spin state of a hydrogen atom is described by its two quantum numbers \FMf), the most general characterization 
would be to write a quantum-mechanical density matrix Pfm f ,f'M' ( r i v )> so that e.g. /9oo,oo(f, v) would describe the number 
of H atoms in the ground state per unit volume, per unit volume in velocity space. The use of the density matrix allows 
for correlations between states of different F or Mf, just as the U and V Stokes parameters allow for correlations between 
the horizontal and vertical polarizations of light in radiative transfer theory. In the realistic high-redshift Universe, we do 
not expect correlations between states with F = and states with F = 1, because their energy splitting implies any such 
correlations are destroyed on a timescale of h/AE ~ 10 -9 s much shorter than the collisional or radiative transition timescale. 
There may however be correlatio ns between states with F = 1 but different Mf, since these are degenerate (except for small 
effects such as Zeeman splitting : ICoorav fc FurlanettdEoolh . 

This general problem can be simplified considerably for the special case where (A) the spin and velocity relaxation times are 
fast compared to the diffusion or free-streaming time across scales where bulk velocity or temperature gradients are significant; 
(B) the spin and velocity relaxation times are fast compared to the expansion, rotation, and shearing timescales of the fluid, 
and the Hubble and Compton-heating timescales; (C) an isotropic radiation field with smooth frequency dependence (such as 
the CMB); and (D) collisional transitions are dominated by spin exchange mechanisms. These are very good approximations 
in most of the Universe (i.e. regions of near mean density) during the dark ages. In these regions, the spin relaxation 
time, considering radiative transitions alone, is T+jAAT-f = 2[40/(l + z)\ kyr, which is much faster than the Hubble time, 
100[40/(1 + z)] 3//2 Myr (note that collisions at high redshift act to decrease the spin relaxation time.) The velocity relaxation 
time, i.e. the timescale for the H atoms to relax to a Maxwellian distribution, is 120 kyr at 1 + z = 40. 3 In regions of near mean 
density, the expansion time is of order the Hubble time, and the rotation and shearing timescales are longer; thus condition 

3 The number given here is the inverse decay constant of the slowest-decaying perturbation from a Maxwellian distribution. Technically, 
it is computed as the inverse of the smallest positive eigenvalue of X^,^ defined in Eq. 1711 . 
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(B) should be valid. Prior to collapse the expectation is that the velocity and temperature are smooth on scales smaller than 
the Jeans length; since the Jeans length is roughly the distance an atom can travel in a Hubble time (neglecting collisions), it 
follows that if (B) is valid, then (A) should be valid in either the diffusion or free-streaming cases. The CMB is isotropic and 
has blackbody frequency dependence, so (C) should be valid unless the 21-cm radiation adds (or subtracts) significantly to 
(or from) the CMB temperature. The me an 21-cm brightness temp erature from high-redshift gas has been computed to be 
< 1 per cent of the CMB temperature fe.g. lLoeb fc Zaldarriagall2004^ . so (C) should be valid in regions of order me an density. 
Final ly, the ratio of dipolar relaxation to spin exchange rate is of order 10~ 2 or less at all relevant temperatures iZygelmanl 
120051) . so (D) is valid. 

The simplifications allowed by (A-D) are as follows. Condition (A) implies that the problem of solving for the density 
matrix and associated master equation can be solved locally, eliminating the position r from the problem. Condition (B) 
implies that we can treat the determination of the distribution function within the steady-state approximation, rather than 
evolving the Boltzmann equation from the recombination epoch forward to the epoch of interest. Condition (C) implies that 
the radiative transition rates are independent of the magnitude and direction of the atom's velocity and are isotropic (i.e. no 
value of Mf is favoured); the smoothness of the frequency dependence is important because otherwise it would be possible 
for atoms moving in different directions to have different 21-cm excitation rates due to the Doppler shift. Condition (D) 
further implies that there is no coupling between the direction of spins and the direction of velocities. Thus in addition to 
being isotropic under rotations of both velocities and spins together, the problem is isotropic under rotations of velocities and 
spins separately. Therefore once a steady-state solution is reached, one can assume that the density matrix Pfm f ,f'm' ( v ) is 
isotropic under rotations of the spins, i.e. it takes the form 

( Mv) \ 

(5) 



Pfm f ,f>m'(v) — 



fi(v)/3 
/i(w)/3 

\ .AH/3 J 

in the { 1 00) , |1 — 1), |10), |11)} basis. Here /f(v) is the level-resolved distribution function, and it represents the number of 
hydrogen atoms in the hyperfine level with total spin F £ {0, 1} per unit volume in physical space, per unit volume in velocity 
space. The total number density of hydrogen atoms is simply 

£ / f F (v)d 3 v. (6) 



nm 

The key problem solved in this paper is to determine the distribution function /f(v) under the above-described conditions 
(A-D) by constructing the Boltzmann equation and finding its steady-state solution, and the associated 21-cm emission. 
Computationally, this is a significant simplification as compared to solving for Pfm f f'm' { r i v )- The limitations of this 
approach should be kept in mind however. The most serious limitation is condition (C), which can be violated if the 21-cm 
line radiation ever becomes comparable in intensity to the CMB (i.e. if the line brightness temperat ure becomes ~ 7%). This 
does not occur in regions of the Universe near mean density, but it may happen in minihaloes (e.g. IShapiro et ai]l2005l) and 
thus our results should be used with caution in this case. 



2.2 Evolution equations 

The distribution function evolves according to several effects including radiative transitions, collisions with various species, 
Lyman-a scattering, and bulk velocity gradients. In general one can decompose the rate of change of the distribution function 
into these effects as 

M«) = ./£ ad) H + /£ o1 - h - h) («) + /|? 0l ' H - He) («) + /i col ' H - e_) («) + /i coliH " H+) («) + .^ LyQ) («) + /^ b) H. (7) 

The term associated with the bulk baryonic fluid velocity Vb is, 

j^w=.-v yy , ( 8 ) 

or 1 ov3 

because the velocity v of the atom relative to the baryonic fluid changes as the atom moves to a region of different Vb even if 
there are no forces acting on the atom. (The Hubble expansion is considered to be part of the velocity gradient.) We neglect 
this term since the associated timescale is of order the Hubble timescale (in mean-density regions) or the collapse timescale 
(in collapsing regions), either of which is long compared to the spin-relaxation timescale. 

Of the remaining terms, the radiative transition term is the simplest to compute. The photon phase space density in the 
21-cm line is given by the blackbody formula M = (e T * / " r "» — l) -1 . We neglect the change in momentum of the atom due to 
absorption or emission of a 21-cm photon, which is an excellent approximation since the recoil energy in temperature units is 
771H = 2 x 10 16 K. Then the radiative term becomes 
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A (rad) H = -/,S rad) W = -A{i+Ar)Mv) + aAMf {v). (9) 

The collisional term is more complicated. For the case of H-H collisions, we have 



i 



/tfoLH-H)^ = K FF ,(v'-v)f F (v)f F ,(v')d 3 v' 

F'=0 

+ ~ Yl J K F , F „(v" -v')f F ,(v')f F „(v")p FlF , F „(v\v',v")d 3 v'd 3 v". (10) 



pt pi, 



Here K F i F n is the product of cross section and relative velocity for two hydrogen atoms in the F' and F" level, and 
Pf\f'f" (v\v' , v") is the probability density for such a collision to produce a final-state H atom in the F level with velocity v. 
The factor of 1/2 in the second term avoids double-counting of collisions since the integral over d 3 v' d 3 v" extends over all of 
R 6 . The probability p F \ F ' F " { v \ v> , v ") must integrate to 2 since the final state has two H atoms: 



E 



Pf\f'F' 



,{v\v' ,v")d 3 v = 2. (11) 



A similar but simpler collision term applies for H-He, H-H + , and H-e~ collisions, but of course in these cases Eq. 1111 integrates 
to 1, and there is no summation over hyperfine levels of the target (i.e. the F' summation in Eq. llUl is not needed). Since He 
is a spin singlet, spin exchange does not occur in H-He collisions, however in general /^ col ' H_Hc ' (w) 7^ because the collision 
changes the velocity of the H atom. For this reason we have included helium, but its effects turn out to be small (< 0.1 per 
cent in the line FWHM and total emissivity). Note that while the direct effect of He collisions is small the 21-cm emissivity 
still depends on the He fraction Yb« because it sets the density of H atoms and the evolution of T^. 

In this paper we will neglect the atom-electron and atom-proton collision terms as we are interested in the period before 
UV or X-ray so urces become impo rtant and begin to produce ionization. There are some free electrons and protons (abundance 
i e ~2x 10~ 4 : Iseager et al1ll999h during the Dark Ages since the cosmological recombination does not run to completion. A 
rough estimate of the importance of the electrons can be obtained by the pre viously computed thermal-averaged spin-change 
(F = 1 — > 0) c ross sections (mov); using the values of the H-H cross section in lZveelmanl I 200,4) and the H-e~ cross section of 
ISmithNl966al) . we find x e {aiov) 1I _ e - /(<tio«}h-h — 0.01 at Tk = 10 K (2 = 21) and less at higher redshifts. We thus conclude 
that electrons change the collisional spin-change r ate at the sub per cent level during the epoch of interest. The spin-change 
rate due to collisions with protons is even smaller JSmitblll966ah . 



2.3 Perturbation theory 

In general the evolution equations are nonlinear and solving them is highly nontrivial. However, in the cases of interest, 
T* = 68.2 mK is much less than either Tk or T 7 and the equations can be linearized by expanding around the thermal 
equilibrium solution. We define 

f { ; h) (v)=n m yF(T k )$ Tk (v). (12) 
to be the thermal distribution function at temperature Tk- 4 We can now define the perturbation 

=/*■(«)-/£%)• (13) 

Since f F h \v) is a constant, all of the evolution equations can be trivially expressed in terms of £f(i>) and its time derivative. 
In terms of these new variables the radiative transition term becomes 
., ^ _*? p - T */ T fc C 1 4. K[\ 4. %hf 

£i rad) W = A 3e - r : /T ^i "hiH.(«) - A(l +SS)£i(v) + 3AN^o(v), (14) 
and working to lowest order in T* this simplifies to 

g*\ v) = \aH±- l) n m <£ Tk (v) - ^AMv) - 3£oH]. (15) 

The equation for the F = levels is related to this by £q (v) = — £1 («)• 

The H-H collisional term in the evolution equation is evaluated by substituting Eq. (11311 into Eq. 110H and expanding 
to leading order in £. In principle, the collision term includes zeroth, first, and second-order terms. However, the zeroth- 
order terms vanish because if the velocities and spin level populations were thermal then the distribution function would be 

4 Note that f F ^ (v) is computed with the spins thermalized at T k . 
© 0000 RAS, MNRAS 000, 000-000 
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independent of time. The remaining terms are 

Cr< H - H V) = -E / K FF ,(v' -v)Z F (v)fW(v')d 3 v' - E f K FF ,{v'-v)f^\v)i F ,{v')d\' 

F'=0 F'=0 

+ E / ^F^K« / '-«oeFK«o/i^ ) («> J P|^(«i« , I t;'')dVdV'. 

F 1 ,F" 

+ E / K F , F ,,{v" - v')£ F/ (v')£ F/ ,(v")p F]F , F „(v\v' ,v")d 3 v' d s v" . 



(16) 



So long as the kinetic and radiation temperatures are 3> T*, we will have |£f(w)| *C f F h \v) and the second-order term on 
the last line of this equation can be neglected. 



3 BASIS FUNCTION METHOD 

It is convenient at this point to expand the perturbation variable in a basis set. We write 

&(v) = E^™<M ?; )> 

n 

where <j> n (v) form a set of basis functions satisfying the orthonormality relation, 

3/2 



m H 



(17) 



(18) 



The normalization has been chosen so that the thermal distribution is one of the basis modes: (j>o(v) = 5>T fc («)■ The remaining 
modes are chosen to be the Gaussian-weighted Hermite polynomials, 



/2a 



(19) 



which indeed satisfy Eq. (1181 . The {<^ n }^L are complete only for spherically symmetric functions (functions that depend only 
on the magnitude v) and in general additional basis functions with angular dependence must be included. However, since our 
problem is spherically symmetric, there is no need to include them. 
Written in terms of basis functions Eq. 1151 becomes 



(20) 



and ^o^ ad ' = — itn^ ■ Evaluation of the collisional rates is aided by defining several integrals. We define for the H-H collisions, 



X 



(col, H-H) 
Fn.F'n' 



(4%a 



2n3/2 



/ E V" (T k )K F , F „ (V" -V)4>n {v)4> {v")& 3 V d V 



+ J y F (T k )K F / F (v - v')<f) n (v)<j} n /(v')(l>o(v)d s v d 3 v' 

- / E yF"{T k )K F , F „(v" - v')p F \ F , FII {v\v' ,v")(t>n,(v)(t) n '(v')(j> {v")d J 'v d^v' dV 

F" 



The first-order terms in Eq. 1161 then can be written compactly as 

/fcol.H-H") irfcol.H-H") 



(21) 



(22) 



Before proceeding, we note that the Maxwellian velocity distribution of Eq. Q assumed in the standard calculation 
corresponds to a perturbation 



Zf{v)\ Ts = {yF(T s )-yF(T k )]n m $T k (v) = -^(-l^Tf 1 - T* l )nm®T h (v), 

or equivalently 
3T, 



16 



(-l) F (T7 1 -T- 1 )n H i<5 rl ,o. 



(23) 



(24) 



In this case only the lowest basis mode (n = 0) is excited. In the full calculation, we will find that in general (p„ can be 
nonzero for any value of n. 
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4 EVALUATION OF H-H COLLISION TERM 

4.1 The relation of K and p to the cross section 

Our objective here is to compute K F 'f"( v " — v') and Pf\f'f"( v \ v ' , v ")- These are related to the differential cross section as 
follows. Defining 

V "+ V ' J // / ^oc\ 

u — and w — v — v , (25) 

we find that 

K F i F n(v" — v') — wa F i F n(w). (26) 
The probability distribution for scattered products is 

p F]F , F „ (v\v\ v") d 3 v = dPF| f F " d 2 f2, (27) 

where dP F \F'F" /dft is the joint probability distribution for the final outgoing angle ft and spin F (which in general depends 
on the angle between ft and w) and 

v = u + -wfl. (28) 

The final relative velocity is w within the approximation of elastic scattering (valid when T* <C Tit) where we neglect the 
hyperfine energy defect in comparison with the kinetic energy. Note once again that the joint probability distribution integrates 
to 2 because there are 2 H atoms: 

±f^^d 2 n = 2. (29) 
p^ J dfl 

The explicit expressions for the differential cross section can be derived by techniques described in e.g. lSmithl l)l966bf) and are 
summarized in Appendix [B] for reference. 

4.2 Evaluation of integrals 

We focus our attention on the third integral appearing in brackets in Eq. 12H : 

T= f J^VF^T^K^p^v'' -v , )pF lF 'F-(v\v\v ,, ) ( fi n (v)<fi n ,(v , )Mv'')d 3 vd 3 v' d 3 v'' . (30) 

F" 

If we can evaluate this integral the first two integrals can also be evaluated by replacing p F \ F ' F " (v\v' , v") with 5 FF / v — v ) 
or 8 FF n8^ (v — v") respectively. In terms of the differential cross sections of i |4.1l this corresponds to setting dPp^p'F" /dfl 
equal to 5 FF i 5^ (ft + w) or S FF iiS^ (ft — w) respectively. 
Upon substituting Eq. 1271 into I, we get 

1= J ^ VF"{Tk)K F , F n{w) ^ F 4>n(u + \ W ^j ^ M ~ ^ W ^j ^° ^ M + \ W ^j ^ U ^ W ^ ft ■ (31) 

Since we are considering only the spherically symmetric modes (fin this integral can be simplified by symmetry considerations. 
We first split the integration over w into integrals over magnitude w and direction w: d 3 w = w 2 dwdw. If we then evaluate 
the d 3 udw d 2 ft integral the result must be independent of the direction of w. We can thus replace w with any unit vector, 
say the unit vector in the third coordinate direction e^. This gives 

J = 47r J ^ y F n (Tk)K F / F /' (w) — ^p 1 ^ (« + </>„, (u - lwe 3 J cfi [u + \we 3 J w 2 d 3 u dwd 2 fl. (32) 

By a similar argument if we decompose ft into spherical coordinates, 

ft — cos 9 e$ + sin 9 (cos <p e.\ + sin ^62), (33) 

then ip can be integrated out to yield 

dP f 
dfl 



J = 8tt 2 y F n (T k )K F > f" (w) — F|F F (fin ~\~ ^(cos9e 3 +sin6>ei)j <j> n , (u - -^ e ^j <Po (u + ^e 3 ) 



F" 

|3„ 



(34) 

Equation 134H contains a 5-dimensional integral over u, w, and 9. In order to solve it efficiently we note that the integrand 
consists of two parts: one is the differential cross section Kf'f"&Pf\f'f" /dfl, which does not depend on u, n, or n'; and the 
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other is the product of basis functions, which does not depend on F, F' , or F" and can be evaluated independently of the 
atomic physics. Therefore we define the integral of three basis functions 



C nn > (w, 6) = I 4> n \u + ^(cos6>e 3 + sin0ei)j <j> n 



W \ , I W \ .3 

u - — e 3 4 u + — e 3 d u. 



r-J-v 2- v- (35) 

It is readily seen that in our choice of basis, Eq. 119L the integrand takes on a special form. The basis function <j> n (v) is a 
Gaussian e~ v l 2a times a polynomial of order 2n in the components (vi,V2,vs) of v. Therefore the integrand in Eq. 1351 is 
equal to a polynomial of order 2(n + n') times the product of Gaussians 



<'xi> \ ~2~S [ u + 2"( cos6,e 3 + sin&ei 



exp 



If w 
'2^{ U ~2 ei 



exp 



2a* \ U+ 2 63 



f 3 r w ~\- 
oc exp < — — ^ |w + — (cos 6ez + sin Oei ) J 



(36) 



where the proportionality constant is independent of u. In 1 dimension the iV-point Gauss-Hermite integration formula gives 
an exact result for products of polynomials of degree < 2N and Gaussians. The same is true in 3 dimensions for the iV 3 -point 
integration formula obtained by performing Gauss-Hermite integration on each of the three coordinate axes. In this case it is 
necessary to choose the weights corresponding to a centroid w(cos Qez + sin (9ei)/6 and la width a/V3, and take N > n + n' . 
The computational demand can be reduced by noting that the integrand in Eq. I|35|l is even under the substitution U2 — > — it2. 
Substituting in Eq. 1351 we obtain 



T=Stv 2 ^2/ F „(T fc )iW'(w)- 



dfl 



-C nn i(w,9) w sin#diud#. 



(37) 



As noted earlier, the remaining two integrals in Eq. 12 II can be evaluated by replacing dPp|F'F" /df2 with Sff/S^^Ci + w) 
and S FF „S (2 \fl - w). This yields 

dPplF/F 



s FF ,s (A> (£i + w) + s FF „s w (n-i 

F" L 

x Cnn' { w , 6) w 2 sin 6 dw d6. 
The S functions of f2 can be converted to functions of 6 according to the standard prescription 
2tvS {2) (n + w) sin6 = S(6-Tv) and 2tt5 (2) (Q, - w) sin6 = 5(8). 
Thus 



df2 



(38) 
(39) 



x (col,H-H) 4 (4 

^Fn.F'n' 1 ^" " 



2\3/2 



5 F f> r^2yF"{T k )K FF „(w)C nn ,(w,n)w 2 dw + / y F (T k )K F , F (w)C nn , (w, 0) w 2 dw 



-2ty 



dP F 



J F n 



dn 



(40) 



Since we are only computing Xp^'p,^ to zeroth order in T* it is permissible to set y F (Tk) — > (2F + l)/4. Making this 
substitution and recalling that K F i F n(w) = wo F i F n(w) we arrive at 



X 



(col,H-H) 
Fn.F'n' 



7r(47T<T 



2x3/2 



6 FF > J Y,( 2F " ' + l)°FF>>(w)C nn ,(w,n)w' A dw + J (2F + l)a F , F (w)C nn ,(w,0) w 3 dw 



- 27 r/^(2F" + l)a F , F „(^ dPF|F ' FA 



dn 



C n „i(w,0)w sin^duidS 



(41) 



5 INCLUSION OF HELIUM 

In reality the primordial baryonic matter did not contain only X H but also 4 He with an abundance n( 4 He)/n( 1 H) « 0.08 
as well as trace quantities (< 10~ 4 ) of 2 H, 3 He, and 7 Li. Collisions with the rare species will be unimportant compared 
to collisions with 1 H but for precision at the level of several percent it is worth considering the role of 1 H- 4 He collisions. 
Helium is an electron spin singlet and hence it cannot undergo spin exchange in collisions with hydrogen; thus in the usual 
approximation of a single spin temperature it does not do anything. If the spin-resolved velocity distribution is nonthermal 
however, collisions with helium can redistribute the velocities of the F = and F = 1 hydrogen atoms and bring them closer 
to the Maxwellian distribution. Indeed if the 4 He/ 1 H ratio were very large, so that collisions with He atoms redistributed 
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the velocities faster than collisions with H atoms could change the hyperfine level populations, then one would expect the 
spin-resolved H velocity distribution to be Maxwellian and the single spin temperature approximation to become exact. 

In this section we expand the Boltzmann equation to include the effects of helium. The first step is to introduce the 
Maxwellian distribution for helium and the appropriate basis modes for describing perturbations of the helium distribution 
function. We may then write the appropriate terms in the Boltzmann equation to account for H-He and He-He collisions. 

In analogy with Eq. 1191 we define helium basis modes 



/ \ n -n— 5/2 -3/2r/ . 1 \|i-l/2 m Hc „ 

(v)=2 n (2n+l)! H 2n +i , 

kBT k v \ ^/fc B r t /m He 



-m He v 2 /2k B T k 



(42) 



which satisfy the orthonormality relation 



0He,n(w)^He,n'( w ) ^ v 



( m Ho 



3/2 



Note that the Maxwellian distribution is simply &nc(v) = ^>He,o(i>), and the thermal helium distribution function is f^' (v) 
^Hci$Hc(u). We may then define the perturbation 



(43) 



(44) 



We now seek a linearized form of the Boltzmann equation for the H-He and He-He terms in analogy to Eq. 141H . The spin 
exchange mechanism does not operate in H-He collisions since He is a spin singlet (S = 0). Therefore the H atom has the 
same value of F before and after the collision, and the H-He collision term in the Boltzmann equation for H is 



;(col,H-He) 

If 



(v) 



J K<-*- H ')(v'-v)f F (v)f H *(v')d a v' + J K^- He \v''-v')f He (v')f F (v'')p^-^(v\v',v'')d 3 v'd 3 v'' 



(45) 



Here 7f' H Ho '(u>) is the product of relative velocity and cross section at relative velocity w, and p^ 1 Hc \v\v' ,v") is the 
probability density for the final velocity of the H atom. Writing this in terms of perturbation variables gives 



i(col.H-He) / n 



K 



(H-He)/ / _ 



(«' - v)[£ F (v)f£>(v') + f^(v)ZMv')}d 3 v 



+ 



K K >[v -v)[f^ e '[v)^ F {v )+£hc(«)/> >{v )]p^ >(v\v,v )dvdv , 



(46) 



where we have used the fact that the zeroth-order terms on the right-hand side vanish. We have also dropped the second-order 
terms in £. 

Integration of Eq. 1146 [I against (j> n (v) gives 



*(col,H-He) 



v /-fcol.H-He) v(col,H-Ho) *(col,H-Hc) 

-nHel2^ Y nn>e Fn > TlHI Fn,He n' £ H e, n ' 



(47) 



where the coefficients are: 



(47T(T 



2\3/2 



K 



(H-He) / / 



{v - v)<j>„ (v)(j> n , (v)(j>He,o(v')d 3 vd 3 v' 



J ^(H Hc)^» _ w ')^ n ( v )^ He>Q ( t ; / )^ n/ ( w ")p^ Uc \v\v',v")d 3 vd 3 v' d 3 



and 



X 



(col, H-He) 
Fn,He n' 



yF{Tk){^-K(T 



2x3/2 



/ 



r,(H-He) / // 



K ( H - H °\v' -v)Mv)Mv)^, n 'iv')d 3 vd 3 v' 



(«)^>He,n' (v')(j> (i>")Ph («!«', V " ) d 3 V d 3 v' d 3 v'' 



(48) 



(49) 



It will turn out that we need only explicitly evaluate Y nn i as ^Fn 'thTn' wm no * an?ec t the hydrogen level populations at 
linear order in perturbation theory because, as we show below, in the steady state £h c n = 0. 

As in our analysis of H-H collisions the most efficient way to evaluate Y nn i is to focus on the second integral in brackets 
in Eq. 14811 . Substitution of p^ ol ' H-Hc ) (v\v', v") — > 5' 3 '(w — v") recovers the first term. The evaluation of the integral is most 
easily accomplished by considering the centre-of-mass and relative velocities rather than v' and v" as variables: 



tuhv +mucV „ , 3 , 3 „ 3 3 

= and w — v — v ; d v d v =d «d 

771H + mft 



(50) 
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The probability distribution for the final velocity v in the scattering event is given by the analogue of Eq. 12711 

(col,H-He)c ] / ll\ ,3 dP ,2A 

pjj ' (v\v ,v )d v = — — d si, 

df2 

and the specific relation between v and CI is now 

v = u H w, 

mn + miic 

which is different from Eq. 1281 since the He and H masses are unequal. The second integral in Eq. (I48H then becomes 



(51) 



(52) 



K Hc) (v" - l/)^ n (T;)</>He,o(w')^n'(y')PH (« I V , v" ) d 3 V dV d 3 "<" 



K 



(H-Hc) 



{w)4> n U + 



nine 



mu + niHc 



CI 0Hc,O U 



mil 



1 m H + rriHe 



m He \ dP 3 3 2A 

-iu — — d tid wd il 



mn + m He / dCl 



m He we 3 \ dP 2 3 2 ~ 
6„/ | u + ; ) — —w d u dw die. (53) 



mn + muc ) dCl 



In the last line we have used spherical symmetry to choose w to be in the third coordinate direction and replace d 3 w with 
4nw 2 dw. If we write CI in spherical coordinates (Eq. 1331 and integrate out the ip direction we then get 

mH B w(cos 9 e 3 + sin 9 e.\ 



u + 



mn + mHe 



<^He,0 « 



mnwe 3 
mu + mn€ 



mu c we 3 
mn + mm 



I 2 = Sir 2 J K^ H -^(w)4> n 

dP 
dCl 

Once again it is possible to separate the purely kinematic terms in this integral from the cross sections. We define 



x ^-r- w 2 sin 9 d 3 u dw d9. 



D nn ,(w,9) = I (f>n 

so that 



u + 



mH e w(cos 9 e 3 + sin 9 ei) 



mu + TTlHc 



!>He,0 U — 



m H we 3 
m H + m Hc 



mn c we 3 \ 3 
u H | d u, 

mn + mm 



1 2 =$tt~ z I K (ll - Ilc) (w)^D nn ,(w,9)w'sm8dwd9. 
dCl 



dP 



(54) 



(55) 



(56) 



Just as for C nn i, it is possible to evaluate D nn i exactly by Gauss-Hermite integration on each of the three coordinate axes. 
However in this case the centroid of integration is at —mnmnaw(cos 9 e 3 + sin 9 ei)/[(mn + mHc)(2mn + mHe)] and the la 
width is ^Jk B T k /(2mn + m Hc ). 

With this value for T2, and recalling that the first integral in Eq. 1481 can be obtained from the second by the replacement 
?4 col ' H-He) W«', v") -> 5 (3) (w - v"), or equivalents dP/dCl -> 6 (2) (Cl - w), we find 



Y n „, = (47T(7 



2n3/2 



8tt / K {H - tic> (w)S^ > (Cl-w)D nn ,(w,9)w 2 sm9dwd9 



-8tt 2 / K (ll ~ Ilc) (w)^-D nn/ (w,e) w 2 sin9dwd9 
J dCl 

Using Eq. 1391 this simplifies further to 



Y nn , = (4n a *) 



2\3/2 



/ K <ll ~~ llc) (w)D nnl (w,0)w 2 dw~2TT [ K (ll -' ac> (w)^-D nn/ (w,9)w 2 sin9dwd9 
J J dCl 



(57) 



(58) 



which is now suitable for numerical evaluation. 



6 SOLUTION OF THE BOLTZMANN EQUATION 

We are now interested in the steady-state solution to Eq. Q. For simplicity we restrict our attention to the case of a purely 
atomic gas of H and He and no Lya radiation. These conditions are expected to pertain to the era before the first sources of 
ultraviolet and X-ray radiation. 

6.1 The solution 

As currently constructed the Boltzmann equation can be written as 

£ X£ S, (59) 
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where X is the relaxation matrix and S is the forcing term. The dimension of the vectors £ and S is 3iV, where N is the 
number of basis modes used. There are N modes each for the perturbations in the hydrogen F = phase space density, for 
hydrogen F = 1, and for helium. The forcing term is purely radiative 



Sin — _ Son — -A ~ ij n HI<5 n ,0, 

and can be read off from Eq. (1201 along with Sne,n = 0. The relaxation matrix is 

Xv(rad) . wfcol.H — H) , w(col,H — He) . wfcol.Hc — He) 

= '+n H iX^ '+n m X '+nmX ■ 



(60) 



(61) 



The various contributions to the relaxation matrix exhibit sparseness patterns when broken into the blocks corresponding to 
H(_F = 0), H(i ? = 1), and He. This turns out to be very useful in solving the Boltzmann equation because only some of the 
rates are relevant. The radiative contribution can be read off from Eq. 1201 and it breaks down as 

/ I -31 \ 

-I 31 , (62) 




r (rad) 







/ 



where I is the N x N identity matrix. The H-H collision term only has non-zero entries for the hydrogen atoms, 



,(col,H-H) 



/ Y (col,H-H) v (col,H-H 
A 00 A 01 



V 



.(col, H-H) 

no 





v (col,H-H) 
A ll 





(63) 



The H-He collision term has a particular structure due to the elastic cross section being the same for H(i ? = 0) and H(F = 1), 
and because the spin-exchange cross section vanishes: 



, (col, H-He) 



/ (nHci/n H i)Y 







,(col,H-He) 
v 0He 
(col, H-He) 



V x 



(n Hc i/^Hi)Y X 1Hc 

(col, H-He) y(col,H-He) y(col,H-He) 

HeO A Hcl HcHc 

(col,H-Hc) 



(64) 



and Y is the matrix of H-He collision integrals defined in Eq. 1481 . 



i wfcol.H — He) \jr 

where X^ c0 ' = X Hel 

Finally, for the He-He term, only the X^"^ 6- block is nonzero. 
The key to solving the Boltzmann equation is to introduce the change of basis from the abundances of H(F = 0) and 
H(F = 1) to the new modes 



^A.n — £l,n — 3^0, n, (,S,n — £l,n + and £ H c,n — ?He,n- 

The new basis {£a,?i> £e,b> £tle n} is related to the old basis {£o,n, £i,n, £He,n} by a matrix, 



€' = R£, R = 



31 I \ 
I I 

) I 



(65) 



(66) 



The relaxation matrices X can be expressed in the new basis via X' = RXR \ and the source is expressed as S' = RS. It is 
at this point that a simplification occurs. Although X is in general non-sparse (cf. Eqs. I62H64H . the first column of blocks of 
X' read as follows: 



Xaa 
Xea 



4^I + ^[3X 



(col, H-H) 
00 



3X 



(col, H-H) 



.(col, H-H) 



v (col.H-H), . v 
X^ ; ] + TCHclY, 



~4~ l 



.(col, H-H) 
V 00 



+ X 



(col, H-H) 
01 



.(col, H-H) 

Ho 



v (col,H-H), 
A ll J 



2 



-lf'-'X^-^, and 



x;, 



= 0. 



(67) 



The latter expression may be evaluated by substituting Eq. I41H in for X^ C °'; H H '. This gives 



-X-En.An' = 7T(47TCT ) — 



(-1) F ' _1 (2F + l)a F , F (w)C nn ,(w, 0) w 3 dw 



-2ty 



E ( _1 ) F _1 ( 2F " + — ,1 F " c nn' (w, 6>) w 3 sin 6> dw d6 



dfi 



(68) 



(The first term in the square brackets in Eq. 14 II does not contribute because it trivially cancels when summed over F' .) Using 
Eq. 1B8L one may verify directly that 
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£ (-lf'-i(2F" + l)a^„H^g^=0 (69) 

for all ?/; and 6 and hence the last integral in Eq. <)68fl vanishes. By integrating Eq. I)69fi over angles, we find 

£ (-l) F ~ 1 (2F" + l)a F , F »H F|FF dft = 0. (70) 

r . /- F,F',F" ''^ 

Noting that the total cross section is symmetrical, ctf'f" (w) = ap"F' (w), and that F" is a dummy variable in the summation 
in Eq. I|70^ . we see that in fact the first integral in Eq. I|68^ also vanishes, so that X En An , — 0. Therefore the matrix X', and 
in particular its first column, reads 

^ X AA X AE X' AHo \ 

X'= X' EE X EHo . (71) 
\ o x HcE X HcHc J 

As we are about to see, we will only need to evaluate the first column of X'. The source vector in the new basis is 
Sa,„ = 3A ( - 1 ) nmSnfi, <S>e,„ = 0, Su c , n = 0. (72) 



Tk 

Thus only one entry of the vector S' is nonzero. 

The particular sparseness patterns of X' and S' axe useful because they imply that the evolution of £ E „ and £H e ,n does 
not depend on £ A: „ an< l nas no sources. Therefore these components of £' - and more generally, the corresponding distribution 
functions fo(v) + fi(v) and fnc(v) - will relax to thermal equilibrium. We have chosen to expand our distribution functions 
around the thermal equilibrium solution at temperature Tj, and hence we will have 

= Ckn = 0. (73) 

This result allows us to simplify the equation for £ An to yield £ A n = — J2 n ' ^"An.An'Ck.n' + ^>A,n- K we consider the steady- 
state solution with n = we conclude that 

& = (Xaa)" 1 ^. (74) 
Numerical evaluation of Eq. I74H can be simplified by exploiting the result X EA = and using Eq. 1671 to write 

X'aa = X AA - X EA = 4|Ul + n H i[X< c ol ' H - H) - X<f' H - H) ] + n Hc iY, (75) 
so that 

&,„ = 3 ^ - l) n m ({4^AI + n H i[X( c ol ' H - H) - Xr'^ H) ] + nu^y 1 ] , (76) 

\ / n0 

where the subscript ()„o denotes the entry in the row corresponding to the nth basis mode and the column corresponding to 
the 0th mode. 



6.2 A subtlety: CMB heating of the gas 

Before continuing, we discuss one subtlety of the steady-state solution to the Boltzmann equation, Eq. The method of 

solution depends on the observation that X EA = 0, which we proved for the elastic approximation cross sections. It is natural 
to ask whether X EA = would still be true if we dropped the elastic approximation and used the full cross sections instead. 
While the answer to this questions is in fact no, and this drastically alters the steady-state solution to Eq. 159H if the full 
(inelastic) cross sections are used, the solution of Eq. (1761 remains the physical solution to Eq. 0. 

In addition to being an approximation for the cross sections, the elastic approximation (in the simple form used here) 
treats the scattering as elastic, i.e. the total kinetic energy of the atoms is the same before and after collision. For this reason, 
the elastic approximation does not include heating of the gas by the CMB in the 21-cm line (the process in which the CMB 
heats the spins, which then collisionally transfer energy to the kinetic degrees of freedom of the gas). Of course, if this process 
is taken into account, the only possible steady-state solution to Eq. 159H is for both the spin and kinetic degrees of freedom 
to come to thermal equilibrium with the CMB. Recall, however, that in the above derivation we dropped the explicitly time- 
dependent terms associated with the Hubble expansion because the Hubble time during this epoch is much longer than the 
radiative or collision times. As we show below, the 21-cm heating time is even longer than the Hubble time and thus, while 
in a nonexpanding universe the spin and kinetic degrees of freedom would have enough time to thermalize with the CMB, in 
our expanding Universe such heating is extremely inefficient. For this reason the solution to the Boltzmann equations derived 
in the previous section by neglecting processes with timescales longer than the Hubble time is the physical solution. 
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A relatively simple argument can be used to estimate the timescale for heating in the 21-cm line. In the spin temperature 
approximation, we know that the net rate of collisional de-excitation per H atom is 

Fde-ox = nm[{(nov)(Tk)ya(T s ) - (a iv)(Tk)yi{T s )] « -nm{<Tiov)(Tk)T i ,(Tj~' 1 - T" 1 ), (77) 

where (<7io«)(7fc) is the thermally averaged cross section for de-excitation, and the approximation makes use of Eq. © and 
the principle of detailed balance (aoiv)(Tk) = 3e _T * ,/Tfc (aiov) {Tk). The heating rate in the 21-cm line is then kBT+Tde-ex per 
atom. Comparing this to the heat capacity for a monatomic gas, C v — 3fcs(l + /h c )/2 per H atom, implies that the heating 
rate is 

• feTJd-, _ nmWiovHn)!*^ 1 -Tr 1 ) f7R , 

Jfe|21cm - p; - , , r , (tO) 

0„ 2(1 + /He J 

implying a timescale for heating of the gas 

t = T k _ 2(1 + fnc)T k . . 

H T fe | 21cm nHi^xovXTfc)^^- 1 -^" 1 )' 

Direct computation with the standard evolution of T s shows that this is always much longer than the Hubble time, e.g. 
th = 50 Gyr versus H' 1 = 100 Myr at 1 + z = 40, and t h = 170 Gyr versus H' 1 = 200 Myr at 1 + z = 25. A more refined 
analysis is possible using kinetic theory but since the heating timescale is so much longer than any other timescale in the 
problem we need not do this calculation. 



7 LINE EMISSIVITY AND PROFILE 
7.1 Determination of emissivity 

Having determined the full distribution function for the various levels of hydrogen, we now turn to the problem of determining 
the 21 centimetre emissivity. This emissivity is most conveniently expressed using the time derivative of the photon phase 
space density Af. The phase space density is related to the familiar specific intensity via Af = c 2 /„ / (2hu 3 ) and is more 
convenient for cosmological calculations because, unlike I v , it is conserved along a trajectory. Its time derivative is 



2hv' A 



(80) 



In this equation j v comes from spontaneous emission, whereas k„ comes from a combination of absorption and stimulated 
emission (the latter giving a negative contribution) . The spontaneous emission term is given by the usual expression 



hvA 
4tt 



dnm(F = 1) 



duii dV 



dun 



hcA 



4tt 



(81) 



where the component of velocity along the line of sight is vu = c(l — v/vio), V is the volume, we have integrated over the 
irrelevant transverse velocity and assumed jl — u/uio\ <§C 1. Since we have calculated the distribution function in the bulk 
rest frame of the gas v\\ is measured in this frame. The stimulated emission contribution is obtained by multiplying this by 
the photon phase space density M and the absorption term is then obtained by replacing /i — > fo and multiplying by 3 to 
account for the lower statistical weight of the F = level: 



~ W 1 " = / LM*0 - 3/o(«)] d 2 v ± . 



We thus find 



87T!/ 3 



Mv)d 2 v ± +Af 



J lfi(v) 



3f (v)]d 2 v ± 



Splitting /f(v) into the thermal piece fp (v) and non-thermal piece £f{v) yields 



Af = 



c 4 A 
Stvu 3 



[yi(T k )+Afyi(T k )-3Afy (T k )]nm / («) d 2 v ± + / ^(v)d 2 v ± +A/ / [^{v) ~ 3£ («)] d 2 v ± 



In the limit where T* T k , T 1 we have 



yi(T k ) +Afyi(T k ) - 3Afy (T k ) = 



(3c 



-T t /T k 



(l + 3e- T */ T fc)(e T */ T T -1) 



Substituting this into Eq. (I84H and re-expressing it in terms of basis modes we find 



Af = 



c A 



£ / _L \ ^ ^ 



(82) 



(83) 



(84) 



(85) 



(86) 
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where ij) n (v\\) = J 4> n (v) d 2 v± is given by the formulas in Appendix fK\ Since £' s n — 0, we have £i„ = and since 

AT « r 7 /T* 2> 1 the third term dominates over the second. Thus we arrive at 



M ■ 



c 4 A 
8-Kvf 



( 1 - 7jr) mnVoO||) + 7jT X/^.nV'n^ll) 



(87) 



The observed brightness temperature Tf, due to the Hi 21-cm signal is obtained by the usual equation, Tf, = hv b s AAf /ks, 
where i/ Q bs = fio/(l + z) is the observed frequency today and A7V is the change in the phase space density. This gives 



hvw 



-AW = 



JVdt = 



t 



A' 



drii 



T 



JVdru 



(88) 



fcjj(l + z) 1 + z j c(l + z)J l + z c(l + z) 2 

where rii is the comoving distance in the radial direction. In this equation A/" should be evaluated using Eq. 1)87^ . Note that 
Eq. 187H depends on vu , which we have defined to be the radial velocity of the H atoms that contribute to the signal, relative 
to the bulk flow of the gas. This velocity is 



H(z) 
l + z 



(89) 



where R\\(z) is the radial comoving distance to redshift z in an unperturbed universe, H(z)/(l + z) is the cosmological velocity 
gradient, and the peculiar velocity ufj pcc ' is needed since the atom's velocity is measured relative to the gas and not to the 
Eulcrian coordinate system. The negative sign is required since Vu is the velocity of the atom relative to the gas (not the other 
way around) although in practice it does not matter since i/>„(i>||) is even. 

7.2 Large-scale perturbations 

For large-scale perturbations where the velocity gradient and density do not vary significantly across the line profile, Eqs. I|880 
and lHUl can be combined to give 

(pcc) l _1 T f 

Afdv t (90) 



T, = 



H(z) dv 
+ 



l + z 



dm 



c(l + z) 2 

If we further use Eq. 18 7H for the line profile, we find 

„(poc)"> - 1 



Th = 



H(z) M 
' + 



3c" AT* 



32nuf (l + zY Um V T:+ 



1 + z ' dr 

where the effective spin temperature Tf 1 is given by 



V(2n + l)!£k,n 
2 n n\ riHi 



(91) 



(92) 



This result should be compared with Eq. ff). Note that for the standard assumption of individual Maxwellian distributions 
for F = and F = 1 (Eq.|24jl, we have &,nk = -3T*(T S _1 - T fc _1 )n,Hitfn,o/4, and hence Tf\ Ta = T s . 

For linear-regime perturbations, one may write Eq. 1911 as a function T b (z, 8 b , dv^ cc ^ /dry ). (We consider T&, Tf 1 , and 
nni to be functions of Si,.) We expand this as 



T b = T b 1 - 



l + z 



d^ pec) 



H(z) dry J + dS b 



(93) 



where T> is the brightness temperature at mean density with no peculiar velocity. Taking the power spectrum of Eq. I|91|l . 
one may obtain tearkana fc Loebll2005h 

2 7-> n \ 4 , 



P Tb (fc) = P M (fc) + (fc) + M 4 P M 4 (fc) 

where n — fc||/fc is the cosine of the angle between the line of sight and the wave vector, 
(fjfij) Ph(k), P^(k)=kT b ^P Sb , Vb (k), and P M 4(fc) =k 2 (T b ) 2 P Vb (k). 



dS b 



(94) 



(95) 



Here Ps b (k) is the power spectrum of the baryon density fluctuations, P Vb (k) is the power spectrum of the baryon velocity 
perturbations, and Ps b ,v b (k) is their cross-spectrum. 



5 iBharadwai fc AlH 120041) and lBarkana fc Loebl l2005fl used the additional relation v b = 6 b /k, which is valid for linear evolution on large 
scales in the matter-dominated era. 
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Figure 1. The temperature shown for mean density intergalactic gas at several redshifts. The solid (blue) curve is the spin temperature 
as a function of atomic velocity v; the long-dashed (red) line is the CMB temperature; the short-dashed (g reen) line is the ki netic 
temperature; and the dotted (black) line is the spin temperature computed by the standard formula (Eq. 23 of lMadau et al.|ll997l with 
y a = 0). Note that the spin temperature is closer to the kinetic temperature for the fastest-moving atoms. 



7.3 Results for high redshift intergalactic gas 

We have computed the 21-cm emissivity and line profile for the high-redshift intergalactic gas. The CMB tem perature was 
deter mined using the usual redshift scaling, T-y = (1 + z)Tq, where Tp = 2.728 K is the present day temperature jFixsen et alJ 
Il99d) . The gas kinetic temperature was determined using Recfast JSeaeer et alJll999l) ; we consider only the era before X-ray 
and shock heating become important in controlling the temperature of the gas. We do not know precisely when this occurred, 
but simulations show t hat even in the abse nce of astrophysical radiation sources, shocked minihaloes dominate the mean 
21-cm signal at z < 20 (|Shapiro et aljEoOflf) . and probably dominate the power spectrum even sooner. We therefore cut off 
our plots at z = 20 but one should be aware that the approximation of adiabatically expanding unshocked gas may break 
down before then. 

In Fig. Q we show the spin temperature of the gas at mean density, as a function of the atom's velocity. For the purposes 
of this figure, we defined a velocity-dependent spin temperature, 

T.{v) = T * , (96) 

ln[3/o(«)//i(v)] 
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which is the generalization of the usual definition T s = T*/ ln(3no/ n i)- We can see that the faster-moving atoms have spin 
temperatures that are closer to the kinetic temperature, which is expected since they undergo more frequent collisions. 
Conversely, the (more typical) slower-moving atoms nearer to the center of the velocity distribution have spin temperatures 
that are closer to the CMB temperature. We thus expect the net effect will be a reduced magnitude of the 21-cm emissivity. 

In Fig. [5] we show how the power spectrum of 21-cm fluctuations on large scales is changed by using kinetic theory in 
place of the usual assumption of a single spin temperature. These results were determined by examining how T b and 8T b /dS 
change with the new calculation and then computing (c.f. Ea . I95[l 

p^{k) _ ( [dT b /d5] \ 2 p^(k) [dT b /os] n i>(fc) _ f % V (Q7 , 

p M o(fc)| Ts \[dT b /dS]\ T J ' p^{k)\ Ts \m/dS\\ Ta %\ Ts ' an i>(fc)k \T b \ T J ■ 1 ' 

The derivatives with respect to the overdensity S were calculated assuming the kinetic temperature varies adiabatically, 
T k oc (1 + 5) 2 ' 3 . As one can see, at the low redshifts z < 60 all contributions to the 21-cm signal (//*, fi 2 , and /i 4 ) are 
suppressed. At higher redshift, the fractional corrections to P M o(fc) and P fi 2(k) become large because dT b /dS crosses through 
zero at z ~ 90. (The zero-crossing occurs because 8T b /dS contains two effects: one is the increase in collision rate and overall 
number of hydrogen atoms in overdense regions, which tends to drive T b to more negative values, while on the other hand the 
increase in Tk means that if collisions are efficient, T„ and T b sho uld go up. The former effect dominates at z < 90, while the 
latter effect dominates at z > 90: see e.g. iBharadwai fc Alil20o4 . 

The 21-cm line profile, <p(v), is needed if one wishes to consider the smallest-scale fluctuations in the 21-cm radiation. 
The line profile is a function of the radial velocity and is proportional to the rate of photon emission, i.e. ip(v) oc N (see 
Eq. 1871 . The line profiles for gas at mean density are shown in Fig. [3] here the line profiles have been normalized such that 
| ^(z) da; | = 1 (where x — v/a). Note that they are wider than Maxwellian, because the fast-moving atoms have a lower 
spin temperature (i.e. farther from T 7 ) than the slow-moving atoms. In Fig. 2] we show the Fourier transform of the line 
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21 -cm line profiles at mean density 

z=55 z=44 




v/a v/a 

Figure 3. The 21-cm line profile shown for mean density intergalactic gas at several redshifts. The solid (blue) curve is the actual 
line profile, whereas the dotted (black) curve shows a Maxwellian profile. Note that the profile is significantly wider than Maxwellian, 
especially at the lower redshifts. 



profile, which determines the Doppler cutoff of the 21-cm power spectrum in the radial direction. In the Maxwellian case, the 
Fourier transform is simply a Gaussian, 

<p(0) ' ( ' 

where kr = (1 + z)~ H(z)a ■ Unsurprisingly, the wider-than-Maxwellian line profile in velocity space manifests itself as a 
cutoff at smaller radial wavenumber ku. This leads to an additional suppression of modes in the 21-cm signal with |feii| greater 
than a few hundred Mpc -1 , which unfortunately makes them even harder to observe. It also means that pro posals to measure 
the I GM temperature using the redshift-space anisotropy caused by the thermal motions of the atoms ijNaoz fc Barkanal 
120051) will have to consider in detail the non-Maxwellian nature of the line profile. The widening of the line profile can also 
be illustrated by plotting the increase in FWHM relative to the Maxwellian case; this is shown in Fig. Qj] The line profile for 
adiabatically compressed gas is shown in Fig. |S| Note that the high-density gas (e.g. in the bottom- right panel with 1 + 5 = 8) 
has a nearly Maxwellian line profile, as expected when collisions dominate over radiative transitions. 
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Fourier transform of line profile 

z=55 z=44 




Figure 4. The Fourier transform of the line profile shown for mean density intergalactic gas at several redshifts. The horizontal axis has 
been converted from the conventional kms -1 to the redshift space wavemimber, kn = (1 + z)~ l H(z)M\\ . 

8 DISCUSSION 

The redshifted 21-cm radiation from the dark ages is a promising future cosmological probe, and in recent years there has 
been substantial progress on understanding the theory of this radiation. However in the analyses to date, it has been assumed 
that the Hi atoms can be described by a single spin temperature T s and a kinetic temperature Tt. We have removed this 
assumption by describing the gas with a full joint spin- velocity distribution function /f(i>)- We find that although the overall 
(spin-summed) velocity distribution is Maxwellian, the individual velocity distributions of H I in the excited and de-excited 
hyperfine levels are not. This leads to two effects: a widening of the 21-cm line, i.e. width greater than \/kBTk/mH', and a 
suppression of the overall 21-cm signal. 

The suppression of the signal on large scales amounts to an effect of up to ~ 5 per cent in the power spectrum. While 
this is small, it will be es sential to consider it if precision cosmological con straints are ever to be obtained from the pre- 
reionization 21-cm signal llLoeb fe Zaldarriagal 12004 : iNaoz fc Barkanal l2005"l. On small scales (~ 10 -2 Mpc) the 21-cm line 
may be broadened up to ~ 1.6 times its Maxwellian width. However, observing these very-small scale features may prove 
to be even more difficult than the large-scale features since the amount of power per mode is less on small scales, where 
P(k) oc fc -3 . An alternative metho d of identifying small-scale structures is to lo ok at 21-cm absorption features in the spectra 
of very high-redshift radio sources dCarilli et alJbooi iFurlanetto fe Loebll2002T) . but it is very unlikely that such sources will 
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70 80 



Figure 5. The fractional change in Vim, the velocity at half-maximum of the line profile, as a function of redshift. Note that in the 
velocity-independent calculation (for a Gaussian line profile) 1^/2 1 t s = V2 In la ~ 1.18<r. 



be found behind a collisionally spin-coupled IGM except in models with unusually early structure formation. Another way 
that the line profile could be important is in regions such as minihaloes that have significant optical depth iShapiro et alJ 
2005), since line self-absorption is suppressed when the line is broadened. In this case the contribution of the minihaloes to 
the large-scale power spectrum could in principle be detectable even if individual minihaloes are not identified, although even 
a small amount of Lya emission could deco uple the IGM spin temperature from the CMB and cause the IGM to overwhelm 
the minihalo signal 1 Furlanetto fc ()l1|20IKii. Note that in the case of a radio source or minihalo, the absorption/emission line 
profile would be modified since one is no longer seeing the gas against a background at brightness temperature T 7 . 

At some point in the evolution of the Universe, UV sources turn on and the IGM temperature is affected by Lya scattering 
as well as by collisional and 21-cm transitions (see iFurlanettolBood for some recent models). A detailed analysis of the spin- 
velocity distribution in this case is beyond the scope of this paper, although we suspect that the single spin temperature 
approximation is very good for this case. This is because, at temperatures greater than a few Kelvin and typical IGM Lya 
optical depths (of order 10 6 ), the colour temperature T c of the Lya photons relaxes to the gas kinetic temperature, and the 
width of feat ures in the Lya spectrum is m uch wider than the Doppler shift frequency ~ VLyuv/c induced by the motions 
of the atoms llChen fc Miralda-Escudll2004) ; thus all atoms see a similar Lya radiation field at similar colour temperature. 
This circumstance would of course lead to a velocity-independent spin temperature if Lya and 21-cm dominate over collisions. 
At temperatures of several Kelvin, obtainable i n low-density regions if reionization starts late, the colour temperature can 
exhibit more complicated behaviour if T s 7^ T fc fairatalbOOfl) . This could potentially produce a non-Maxwellian line profile 
and affect the total emissivity; in order to solve this problem one would need to simultaneously solve the radiative transfer 
and Boltzmann equations for the UV photon spectrum and Hi velocity distribution. However the single spin temperature 
approximation would apply to high accuracy if X-ray heating drives Tk to tens of K or more, or if the Lya spin-flip scattering 
rate ever becomes much larger than the 21-cm transition rate (which would thermalize all spins at 2*.). Collisions with electrons 
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21-cm line profiles at z=33 for various overdensities 

z=33; 5=0 z=33; 6=1 
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Figure 6. The 21-cm line profile shown for intergalactic gas at z = 33. We show values ranging from mean density (<5 = 0, left panel) 
through 8 times mean density (<5 = 7, right panel), assuming a kinetic temperature from adiabatic evolution, Tj. oc (1 + <5) 2 / 3 . The solid 
(blue) curve is the actual line profile, whereas the dotted (black) curve shows a Maxwellian profile. Note that the line profile is more 
nearly Gaussian at high densities and temperatures where collisions are more effective. 



may also be important if there is partial ionization due to X-rays (|Kuhlen et alJ .2006): this hyperfine transition mechanism 
is expected to be independent of the atom's velocity because this is negligible compared to the velocity of the electron. 

We conclude that the traditional assumption of a single spin temperature describing the high-redshift gas, while a good 
first approximation, is incorrect at the several per cent level for the mean signal and at the several tens of per cent level for 
the line profile. Any calculation aiming for higher precision must thus include the full kinetic theory analysis and the joint 
spin-velocity distribution of the neutral hydrogen atoms derived here. 



ACKNOWLEDGMENTS 

We thank S. Furlanetto for insightful comments on a draft of this paper. CH is supported in part by NSF PHY-0503584 
and by a grant-in-aid from the W. M. Keck Foundation. KS is supported by NASA through Hubble Fellowship grant HST- 
HF-01191.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for 
Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. 



© 0000 RAS, MNRAS 000, 000-000 



Spin-resolved atomic velocity distribution 21 



REFERENCES 

Allison A. C, Dalgarno A., 1969, ApJ, 158, 423 
Babich D., Loeb A., 2005, ApJ, 635, 1 
Barkana R., Loeb A., 2005, ApJ, 624, L65 

Bharadwaj S., Ali S., 2004, MNRAS, 352, 142 

Bowman J. D., Morales M. F., Hewitt J. N., 2005, preprint ( astro-ph/0512262 ) 

Carilli C. L., Gnedin N. Y., Owen F., 2002, ApJ, 577, 22 

Chen X., Miralda-Escude J., 2004, ApJ, 602, 1 

Chuzhoy L., Shapiro P. R., 2005, preprint (astro-ph/0512206 ) 

Cooray A., Furlanetto S., 2005, MNRAS, 359, L47 

Dalgarno A., 1961, Proc. R. Soc. London A, 262, 132 

Edmonds A., 1960, Angular Momentum in Quantum Mechanics (Princeton: Princeton University Press) 
Field G., 1958, Proc. I. R. E., 46, 240 
Field G., 1959, ApJ, 129, 536 

Fixsen D. J., Cheng E. S., Gales J. M., Mather J. C, Shafer R. A., Wright E. L., 1996, ApJ, 473, 576 
Frye D., Lie G. C, Clementi E., 1989, J. Chem. Phys., 91, 2366 
Furlanetto S. R., Loeb A., 2002, ApJ, 579, 1 
Furlanetto S., 2006, preprint (astro-ph/0604040 ) 
Furlanetto S., Oh S. P., 2006, preprint ( astro-ph /0604080 I 

Gradshteyn I. S., Ryzhik I. M., 1994, Table of integrals, series and products, 5th. ed. (New York: Academic Press) 
Hirata C. M., 2006, MNRAS, 367, 259 

Iliev I. T., Shapiro P. R., Ferrara A., Martel H., 2002, ApJ, 572, L123 

Iliev I. T., Scannapieco E., Martel H., Shapiro P. R., 2003, MNRAS, 341, 81 

Jamieson M. J., Dalgarno A., Yukich J. N., 1992, Phys. Rev. A, 46, 6956 

Jamieson M. J., Dalgarno A., Zygelman B., Krstic, P. S., Schultz D. R., 2000, Phys. Rev. A, 61, 014701 
Kolos W., 1967, Int. J. Quantum Chem., 1, 169 

Kolos W., Szalewicz K., Monkhorst H. J., 1986, J. Chem. Phys., 84, 3278 

Kuhlen M., Madau P., 2005, MNRAS, 363, 1069 

Kuhlen M., Madau P., Montgomery R., 2006, ApJ, 637, LI 

Loeb A., Zaldarriaga M., 2004, Phys. Rev. Lett., 92, 211301 

Madau P., Meiksin A., Rees M., 1997, ApJ, 475, 429 

McQuinn M., Zahn O., Zaldarriaga M., Hernquist L., Furlanetto S. R., 2005, preprint ( astro-ph/0512263 ) 

Morales M. F., Bowman J. D., Hewitt J. N., 2005, preprint ( astro-ph/0510027 ) 

Naoz S., Barkana R., 2005, MNRAS, 362, 1047 

Pritchard J. R., Furlanetto S. R., 2006, MNRAS, 367, 1057 

Santos M. G., Cooray A., Knox L., 2005, ApJ, 625, 575 

Seager S., Sasselov D. D., Scott D., 1999, ApJL, 523, LI 

Shapiro P. R., Ahn K., Alvarez M. A., Iliev I. T., Martel H., Ryu D., 2005, preprint (astro-ph/0512516 ) 

Sigurdson K., Furlanetto S. R., 2005, preprint (astro-ph/0505173 ) 

Smith F. J., PI. & Space Sci., 1966a, 14, 929 

Smith F. J., PI. & Space Sci., 1966b, 14, 937 

Tang K. T., Yang X. D., 1990, Phys. Rev. A, 42, 311 

Wang X., Tegmark M., Santos M., Knox L., 2005, preprint ( |astro-p h/0501081 ) 
Wouthuysen S., 1952, Astron. J., 57, 31 

Zaldarriaga M., Furlanetto S. R., Hernquist L., 2004, ApJ, 608, 622 

Zygelman B., Dalgarno A., Jamieson M. J., Stancil P. C, Phys. Rev. A, 67, 042715 

Zygelman B., 2005, ApJ, 622, 1356 



APPENDIX A: PROPERTIES OF THE BASIS FUNCTIONS 

Here we summarize some useful properties of the basis functions defined in Eq. 1191 . These are needed in order to compute 
the two quantities of observational interest, namely the velocity-marginalized Hi level populations {no,ni} and the 21-cm 
line profile (which in general can be non-Gaussian). 

Both the level population integrals and the line profile are related to the 3-dimensional Fourier transform of <f> n , namely 
<f>n(M) = f (f> n (v)e~ lM d 3 v. This quantity can be evaluated by switching to spherical coordinates for v: 
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4> n (M ) = / 4> n (v)e 1 v cos 9 v 2 sin 9 dv d6 dip = Aty I <j) n (v) "*" v "" — V 1 dv. 
' ' Mv 



i(Mv) 



Recalling that a — ^JkBTk/mn, we then have, using Eq. Ipjl 
2 -™- 1 /2 7r -i/2 [ ( 2n + !),]-!/ 



6„(M) = 



a 2 M 



, V\ _„2 /2<r 2 
#2n+l ( -J e ' 



sin(Mii) dv. 



(Al) 



(A2) 



We note that H 2n +i is odd so the integrand here is even; thus we may replace J °° — > | f^,- Also sin(M?;) = Sse lM " where 
9 denotes imaginary part. If we further substitute x — v/a, we get 



2 -u-3/2 7r -l/2| (2n+ 1) , ] -l/2^ yoo 

2-"~ 1 [(2n + l)!]~ 1/2 
aM 



Q / _H" 2 „+i(x)e 



2 II iMo 



dx 



(A3) 



where in the s econd equality we have used the Fourier transform of a product of a Hermite polynomial and a Gaussian (e.g. 
Eq. 7.376.1 of lGradshtev n fc Rv zhik 1994). Using the specific equation for H 2n +i, we can see that the integrated line profile 



MO) = (-l)"2- n - 1 [(2n + l)!]- 1 / 2 lira ^±±Ml = (-i)^""' x [(2n + l)!]- 1/2 i/ 2lI+ i(0) 

M— >o alvl 

according to l'Hopital's rule. The derivative H 2n+1 (0) can be determined from the series expansion of H2n+i and gives 



(A4) 



, ir (2n + 2)! - V(2n + 1)! 

H 2n+1 (0)-(-l) (n + 1) , -» 0„(O) = 2 „ n , . (A5) 

In addition to the integrals of the basis modes, provided by Eq. 1A5L we also need the contribution of each basis mode 
to the line profile, defined by 



i>n(v\\) = J i(»)d 2 t)i, 



(A6) 



where v« is the component of v along the line of sight and v± is the component in the plane of the sky. The Fourier transform 
of the line profile is trivially determined, 



M M \\) 



V>n(«||)e iA Vn duy = / ^(wK^IHI dun d 2 v± = MM\\) 
2-™- 1 [(2n + l)!l- 1 / 2 



-(_l)» e - CT M II /2 H 2 „ +1 (<tM,|). 



(A7) 



It is in fact this Fourier transform t[! n {M«) that we will need in order to obtain linear theory power spectra. However it is 
conceptually useful to plot the actual function -i/) TO (t)||), so we calculate this here. It is 

(A8) 



P OO P OO P oo 

ipn(v\\) = / <j>n ( J vf, + v^) 2nv ± dv ± = 2n / <p n (v) v dv = 2n / cj> n (v)vd 
Jo v v ii / V| i j 



Here the second equality involves the change of variables from v± to v = ^Jv 2 + Vj_, and the third equality holds because the 
integrand is odd and hence the integral from vu to contributes nothing. Plugging in the specific form for the basis modes 
from Eq. gives 

1 f°° „ fV\ _„2 /2 2 , 1 



1pn(v\\) 



2^/2ir(2n + l)\a 2 



fv\ _..2/,_a 

tl 2n + l y— I f 



v '^~dv = 



2^2^(271 + 1)! a J v ,,/ 



H 2n+ i(x)e x 12 dx 



(A9) 



where we have again substit uted x = via. The integral may be solved by using the Hermite polynomial recurrence relations 
(Eqs. 8.952.1 and 8.952.2 of lOradshtevn fc Rvzhik1ll994h in the form 



Ax 



2 2 i~ 1 j\ 



— H 2j (x)e 



-x 2 / 2 



H 2 j-i{x) H 2j+1 {x) 



-x 2 - /I 



22j-ij! 



2jH 2j - 1 (x)- -H 2j+1 (x) 



-x 2 / 2 



-x 2 /2 



2 2 <J- 1 1 (j - 1)! 2 2 ij\ J " ^ A1 °^ 

for j > 1. In the special case of j = 0, it is trivial to see that the same expression is valid if one formally defines H-i(x)/[(— 1)!] = 
0. Summing this equation from j — to n, we then get 



d.r 



1 

E 237^1-7 #2, (*)e- 



x 2 / 2 



J=0 



H 2n+1 (x) _ x 2 /2 
2 2n n\ 6 



(All) 
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With this equation, we can now solve the integral in Eq. 1A9L 

n n . 

^ +1 (,)e-^ 2 / 2 d^-2 2 "n!^-^^(,)e- 2 / 2 r_ = £ ^H 2j 3) e~^' , (A12) 



which implies 

Wt»(«ii) = , , > r-^ — e II' . (A13) 



This formula is suitable for numerical evaluation and has been used to generate our plots of the 21-cm line profile. 

APPENDIX B: CROSS SECTIONS 
Derivation of H-H Cross Section 

We determine the H-H scattering cross section according to the elastic approximation jDalgarnolll96lh . In this approximation 
we ignore the hyperfine energy defect and treat the collision as a scattering problem with separate potentials Vo(R) and Vi(R) 
for the electron spin-singlet (S = 0) and spin-triplet (S = 1) states. The difference between these potentials results in a change 
in the electronic spin of a particular hydrogen atom, and hen ce a change in its total (nuclear plus electronic) spin angular 
momentum. The derivation below is based on th e formalism oflSmithl Jl966al lbh. although we use the \SIFtM Ft ) basis rather 
than diagonalizing in S and Ms as was done bv lSmithl Jl966ah . 

There have been many other previous determinations of the cross section for scattering of two hydrogen atoms, some 
going beyond the elastic approximation; however the published r esults do not provide the f ull spin and angular dependence 
of th e cross section, tabulating instead the spi n-flip cross section jA llison fc Dalgarnolllflfiflt IZvgelman et alJl200.4 Izvgelmanl 
120051) or moments of the angular distribution lljamieson et alJ|l992l . I2OO0I) . 

We are interested in the differential cross section for hydrogen atoms in the F' and F" hyperfine levels to scatter and 
leave a hydrogen atom in the F level moving in direction fi. It will be assumed that initially the F" atom is moving with 
velocity wez relative to the F' atom, so that 9 = arccos is the scattering angle. In our problem the hydrogen spins are 
unpolarized because the radiation field is isotropic and because in the elastic approximation the atoms cannot be polarized 
by collisions. Therefore we consider only the cross section summed over final magnetic quantum number M F and averaged 
over initial M F i and M F n . 

We denote the total spin by Ft (the vector sum of F' and F"), and use the labels a and 6 for the nuclei. The incident 
wave function of the two hydrogen atoms is then 

j* in ) = ApA e (\ls a ls b }\F'F"F t M Ft )e lkz ) , (Bl) 



where (X, Y, Z) — Rb — R a is the relative position vector of the two nuclei, k = mnw/2h is the wavevector for reduced mass 
7tih/2, I Isq, 1st) is the orbital wave function of the two electrons (with ls a and lsb representing the Is orbitals associated with 
nuclei a and b), \F' F" FtM Ft ) is the spin state of the nuclei a and b and electrons 1 and 2, and Ap $e are the antisymmetrization 
operators for the protons and electrons. We have chosen our coordinate system such that the relative velocity is along the 
third coordinate axis. Also the notation \F a \Fi > 2FtM Ft ) refers to the angular momenta F a i = I a + Si and Fta = It + Sz- 

The antisymmetrization operators in Eq. 1BH commute with the total electronic and nuclear angular momenta S and /, 
but not with F a \ and Fb2- It is thus convenient to re- write Eq. ijBl^ in the basis of eigenstates of S and I (which we denote 
\SI,F t M Ft )): 

|* in ) = ApAe (^\ls a lsb)Y^{SI,FtM Ft \F'F''F t M Ft )\SI,F t M Ft )e' kz ^j 

= I^[|l Sa lst,) + (-l) s |l S asa)](S7,FtM J? jF'^' , i^M Ft >|SI ) F t Mi, t >[e^ + (-l) s+/ e- i ^]. (B2) 
s,i 

The scattered wave function is related to the ingoing wave function by the scattering amplitudes f(fl). These are defined 
in terms of the large-/? asymptotic form for the wave function, 

ikR 

i*)~ e!fcZ + V /( " } ' (B3) 

in which an incident state |$i n ) ~ e !fcz is mapped into an outgoing scattered state |^ SC at) ~ R~ 1 e lkR f(£l). In our case, the 
ingoing wave function has both singlet and triplet parts with different potentials and hence different scattering amplitudes. 
Thus we must superpose singlet and triplet solutions in order to obtain the wave function corresponding to Eq. 1B21 . Since 
our particles are identical, we must also superpose both the ingoing wave function with e lhz dependence, and one with e~ lhz 
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dependence which can be obtained from Eq. 1B3I by a parity transformation. Thus the wave function of Eq. 1B2I yields an 
outgoing scattered state 



1 1' scat) 



27? 



Y,[\^als b ) + {-l) s \ls b ls a )]{SI,F t M Ft \F'F"F t M Ft )\SI,F t M Ft )[f s (n) 



\S+I 



M-n)], 



(B4) 



where the amplitudes /o(f2) and fi(fl) pertain to electron spin singlets and triplets respectively. These scattering amplitudes 
are determined by solving the Schrodinger equation by the usual partial wave method and we describe the potentials used 
and other pertinent details at the end of this Appendix. 

The differential scattering cross section to put nucleus a and electron 1 into an H atom with total angular momentum 
F in the final state is then obtained by projecting l^scat) onto all outgoing states with F a i = F and arbitrary Fb2, and then 
summing the square norms: 



da 
dtl 



(F al — F) = R 2 Y,\ ( F Fb2F t M Ft \^t) \ 



is 



^2(SI, F t M Ft \F' F" F t M Ft ){FF b2 F t M Ft \SI, F t M Ft )[f s {tl) + {-1) S+I f s (-fl) 
s,i 



(B5) 



(In principle we must also sum over the final values of the Ft and M Ft quantum numbers; however since we have neglected 
interactions involving the spin these will be conserved and we may simply use their initial values.) To get our final result for 
the differential cross section for producing a hydrogen atom with total spin F moving in the direction f2, we must average 
over the initial values of the quantum numbers Ft and M Ft in their statistical ratios, and then multiply by 4 because the 
final-state hydrogen atom could contain either nucleus (a or b) and either electron (1 or 2). The summation over M Ft is trivial 
since none of the inner products or amplitudes inside the sum depend on it: 



dP F 



a pi pi, 



dfl 



{2F> + 1)(2F" + 1) 



x[/ s (f2) + (-l) s+ 7 s (-n)] 



E ( 2Ft + V > 



^<S7, F t M Ft \F'F"FtM Ft )(FF b2 F t M Ft \SF F t M Ft ) 
s,i 



(B6) 



The evaluation of Eq. (IB6t is a straightforward but tedious exercise in angular momentum recoupling theory. The 
recoupling coefficients {SI , F t M Ft \F' F" FtM Ft ) are given in terms of the 9-j symbol by 



1/2 1/2 S 

{SI, F t M Ft \F'F"F t Mp t ) = V(2S + 1)(2J + 1)(2F' + 1)(2F" + 1) { 1/2 1/2 / 

F' F" F t 

jEdmondslll96oh . Using the specific values of the 9-j symbols, we find 



o"oo- 

O"10" 

I 

CTll- 
I 

Tll- 



dPou 



dfl 

dPi|oo 
dfl 

dPoioi 

dfl 

dFiioi 

dfl 

dPono 

dn 

dPi|io 

dfl 

dPom 

dfl 

dPi|u 

dfl 



^|/oo(n)+3/ii(f2)| 2 , 
16 

A|/„o(f2)-/n(f2)| 2 , 

^|/ ()1 (n) + /io(n) + 2/ 11 (ri)| 2 , 

^|/or(n) +/io(0) - 2/n(f2)| 2 + i|/oi(fi) - fia(to)\ 2 , 

^|/ ()1 (n) + / 10 (n)-2/ 11 (n)| 2 , 

i|/ i(f2) +/i (n) +2/u(f2)| 2 + g|/ i(O) - /io(f2)| 2 , 



^l/ 00 (n) - /n(n)| 2 + ^|/oi(n) - / 10 (n) 



and 



^13/ooCn) + /n(o)| 2 + ^|/oi(n) - /io(n)| 2 + ^|/oi(n) + /io(fi)| 2 + §|/n(A) 



(B7) 



(B8) 



where f S i(Sl) = f s ((l) + {-l) b+1 /s(-H 
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Derivation of He-H Cross Section 

The He-H cross section is actually simpler to derive than H-H because the He atom has no electronic spin or orbital angular 
momentum. Thus there is only one relevant electronic state, X 2 E + , and there is no spin exchange. The scattering cross section 
is independent of spin and given by the usual formula, du/dA = /(A)| 2 , where /(A) is the scattering amplitude in the X 2 E + 
potential. 



Potentials, Phase Shifts and Scattering Amplitudes 

In order to calculate the cross sections as outlined above we must calculate the scattering amplitude f(fl) in the atomic 
potential V(R). To do this we solve the Schrodinger equation 

-— V 2 * + (V-£)* = 0, (B9) 

where ^S> = (R\^), n is the reduced mass of the two colliding atoms, and E — h 2 k 2 /2/i for an unbound state with wavevector 
k. For a spherically symmetric potential the wavefunction can be expanded in a basis of angular momentum eigenfunctions 

as 

oo 

A) = Y, ^N(R)P N (k ■ A) (BIO) 

where N is the orbital angular momentum quantum number, the ajv are constant coefficients, Pn are Legendre polynomials, 
and the partial wave i/jn satisfies the radial equation 

d 2 ^ N , r, 2 N(N + 1) 2ft. 



lp N =0. (Bll) 



dR 2 

For R>TZ, where V(TZ) < E - N(N + l)h 2 /(2fiTZ 2 ), ip N takes the approximate form (exact for R -> oo or V -> 0) 

ipM — 5jv(fc_R) cos Sn + CN(kR) sin <5jv (B12) 

where <5jv is the phase shift of the iV th partial wave, and SN(kR) and Cjv(fcJ?) are Riccati-Bessel functions, which are related 
to the spherical Bessel and spherical Neumann functions by Sn{x) = xjn{x) and Cn(x) — —xun(x). Using Eq. l|B12fl and 
the asymptotic forms 

lim S N (kR) -v sin [ kR — and lim C N (kR) -> cos [ kR — , (B13) 



it is straightforward to prove the exact integral identities 
2/t 



sin5jv = --^/ V{R)iP N (R)S N (kR)dR, (B14) 



(i 



and 

cosSn = H| [ ( V (R) + ^ I ip N (R)S N -i(kR) dR. (B15) 



h 2 J \ 1 ' R 2 

We use these integral forms to solve for 8n = <5jv(fc) while simultaneously solving Eq. IBllll for tPn(R). Given these phase 
shifts the scattering amplitude, 

1 oo 

/(A) = - V (2JV + l)e i4jv sin<5jv Pjv(fc ■ A) , (B16) 
k ' 

jv=o 

can immediately be found. 

For this calculation we needed the scattering amplitude in the singlet and tri plet b 3 T,t potentials for H-H sca ttering, 



and in the He-H X 2 T, + scattering potential. We used the same H-H potentials used in lSigurdson fc Fmdjm e1 t t(jj 200.4). which 
is a combination of several previously publis hed H2 surfaces, interpolations, and asymptotic formulas jKoloJl967tlKolos et alJ 



is a combination 01 several previously published M2 surfaces, interpolations, and asymptotic lormulas (Kolos 19b/; Kolos et 
Il986t l iFrve et al.lll989lljamieson et al)ll992ll . For the He-H potential, we used the fitting formula bv lTan^^^^na(^99Cjl) 
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